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Abstract 

Multi-target tracking is mainly challenged by the nonlinearity present in the measurement equation, 



introduces a grid-based model in which the state captures target signal strengths on a known spatial grid 



o 

(N 

^ • and the difficulty in fast and accurate data association. To overcome these challenges, the present paper 

, (TSSG). This model leads to linear state and measurement equations, which bypass data association 

and can afford state estimation via sparsity-aware Kalman filtering (KF). Leveraging the grid-induced 

'■ sparsity of the novel model, two types of sparsity-cognizant TSSG-KF trackers are developed: one 

' effects sparsity through ^i-norm regularization, and the other invokes sparsity as an extra measurement. 

^ I Iterative extended KF and Gauss-Newton algorithms are developed for reduced-complexity tracking, 

along with accurate error covariance updates for assessing performance of the resultant sparsity-aware 

, state estimators. Based on TSSG state estimates, more informative target position and track estimates 

OO ' can be obtained in a follow-up step, ensuring that track association and position estimation errors do not 

00 , 

' propagate back into TSSG state estimates. The novel TSSG trackers do not require knowing the number 

• , of targets or their signal strengths, and exhibit considerably lower complexity than the benchmark hidden 

■ Markov model filter, especially for a large number of targets. Numerical simulations demonstrate that 

sparsity-cognizant trackers enjoy improved root mean-square error performance at reduced complexity 
when compared to their sparsity-agnostic counterparts. 
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I. Introduction 

Target tracking research and development are of major importance and continuously expanding interest 
to a gamut of traditional and emerging applications, which include radar- and sonar-based systems, surveil- 
lance and habitat monitoring using distributed wireless sensors, collision avoidance modules envisioned 
for modem transportation systems, and mobile robot localization and navigation in static and dynamically 
changing environments, to name a few; see e.g., ||4l, [Ol, and references therein. 

At the core of long-standing research issues even for single -target tracking applications is the nonlinear 
dependence of the measurements on the desired state estimates, which challenges the performance of 
linearized Kalman filter (KF) trackers, including the extended (E)KF, the unscented (U)KF, and their 
iterative variants H, ||9]- This has motivated the development of particle filters (PFs), which can cope 
with nonlinearities but tend to incur prohibitively high complexity in many critical applications. For multi- 
target tracking, data association has been another formidable challenge, especially when the ambient 
environment is cluttered, and the sensors deployed are unreliable. This challenge amounts to determining 
the target associated with each measurement, where the noisy measurements typically reflect the candidate 
target locations acquired through signal detection in gated validation regions; see e.g., ||3l, Q. Once 
data association is established, targets can be tracked separately using the associated measurements, in 
conjunction with track fusion for improved accuracy. 

The present paper investigates the multi-target tracking problem whereby the available measurements 
comprise the superposition of received target signal strengths of all targets in the sensor field of view. 
Sensors collecting these measurements are not necessarily radars or high-cost receivers, but can be general- 
purpose radio units employing simple energy detectors. The measurements are nonlinearly related to target 
locations but no data association issues arise, because conventional range-gate operations have not yet 
been employed to detect, separate, and localize the targets of interest ||3l. To cope with the nonlinearity 
issue, this paper introduces a grid-based dynamical state- space model in which the state describes signal 
strengths of targets traversing a preselected spatial grid (TSSG) of the tracking field. Because the locations 
of grid points are preset and known, both the measurement and state equations become linear. Further, 
data association is avoided by dynamically tracking the TSSG values rather than directly producing the 
target tracks. Based on TSSG tracking however, data association and track trajectory estimation can be 
performed as a follow-up step, whereby track association and estimation errors do not propagate back to 
the TSSG tracker. 

Similar ideas on bypassing data association at the price of tracking "less informative" estimates have 
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been exploited in recent multi-target tracking schemes, such as the probability hypothesis density (PHD) 
filter II22I . ll30l and the Bayesian occupancy filter (BOF) |[T4l . The PHD filter tracks the so-termed 
target intensity, while the BOF tracks the probability of a grid point being occupied by any target. A 
main advantage of the grid-based TSSG tracker here is that state estimation becomes possible via KF 
applied to linear state and measurement models, at considerably reduced computational burden relative 
to the complexity incurred by the PHD and BOF. Further, the TSSG tracker is novel in exploiting the 
sparsity present in the grid-based state vector, which allows one to leverage efficient solvers of (weighted) 
least-squares (LS) minimization problems regularized by the £i-norm of the desired state estimate. 

Sparsity-aware estimators have been studied for variable selection in static linear regression prob- 
lems, and have recently gained popularity in signal processing and various other fields in the context 
of compressive sampling (CS); see e.g., ||5l, lITTI . 11211 . However, few results pertain to the dynamic 
scenario encountered with target tracking. When measurements arrive sequentially in time, a sparsity- 
aware recursive least-squares scheme was reported in lH], but its tracking capability is confined only to 
slow model variations; see also ||2| for a sparsity-cognizant smoothing scheme which nevertheless does 
not lend itself to filtering; as well as |[29l . where a so-called KF-CS-residual scheme is reported for 
tracking slowly varying sparsity patterns. Different from existing alternatives, the present work develops 
sparsity-aware trackers along with their error covariances, without requiring knowledge on the number 
of (possibly fast-moving) targets or their signal strengths. 

The rest of the paper is organized as follows. Section |II] develops the novel grid-based sparse model, 
for which a sparsity-agnostic KF tracker is introduced in Section |llll Two sparsity-cognizant trackers are 
presented in Sections JV] and |Vl Target position estimation and track formation is detailed in Section |Vl] 
Numerical results are presented in Section IVIII followed by concluding remarks in Section IVIIII 

n. Grid-based State Space Model 

Consider the problem of tracking M moving targets using active (e.g., radar) or passive (e.g., 
acoustic) sensors deployed to provide situational awareness over a geographical area. Targets emit power 
either because they passively reflect the energy of other transmitters such as radar, or, because they are 
active sources such as cell-phones or transmitters mounted on smart cars. Associated with each target, 
say the mth one, is its position vector p^™^ per time k, and the signal of strength s*^'") that the target 
reflects or emits. Sensor n measures the superposition of received target signal strengths, 

M 

2/n,A:= ^/l(C^")s("'^+Z^n,fc, n=l,...,N, k = 1,2, . . . (1) 
m=l 
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where /i(-) denotes the distance-dependent propagation function; d^"^" := ||p[!"'' — qnib is the distance 
between the known position q„ of sensor n and the unknown position vector p[™^ of target m; and Un^k 
is zero-mean Gaussian noise at sensor n. Function h{-) satisfies /i(0) = 1, is non-negative, decreasing, 
and is either assumed known from the physics of propagation or acquired through training [llTI . 

At each time k, a centralized processor has available the measurement vector := [yi^k^ • • • > yN,k]^^ 
based on which the target positions {Pk^^}m=i t)e tracked. Note that the measurement model ([1} 

differs from the one typically considered in radar applications, where a measurement either comes from a 
single target or a clutter, usually in the form of position information obtained from range gate operations 
Q. Each measurement in ([1} comes from a sensor, and comprises the superposition of received signal 
strengths emitted by or reflected from all targets in the sensor field of view. This model considers the 
localization and tracking problems jointly, and avoids the measurement-target association issue. 

One major challenge in tracking and locahzation problems is that the measurements in ([T} are nonlinear 
functions of the wanted target position vectors. A neat approach to arrive at a linear measurement model 
is to adopt a set of G (possibly regularly spaced) grid points at known positions {gj}iLi> where target(s) 
could be potentially located; see also e.g., |[T4l . lITTI . and Q. Using a sufficiently dense grid, it is possible 
to capture the target locations at a prescribed spatial resolution using a G x 1 vector having most 
entries equal to zero except for the {i^'"^}m=i entries given by {a;^*'' ^}m=i' which represent the target 
signal strengths at time k if and only if the m-th target is located at the 4"^^"'^'^ S^^*^ point, that is 
Pfc = Note that if target m is located exactly on a grid point i^. , then x^*" = s^'^' / will 

be the only nonzero entry of corresponding to this target. However, to account for target presence off 
the preselected grid points, it will be allowed for the unknown target signal strength s^*") to "spill over" 
grid points around i^™"* and thus render nonzero a few neighboring entries of x^. Let Q^'^ denote the 
spill-over region on the grid corresponding to target m at time k, such that 7^ is associated with 
s^"^\ \/i G Gk^^- The following assumption on this target occupancy model is imposed: 
(asl) Each grid point i can be occupied by at most one target m at any given time k. 
This assumption can be easily satisfied in practice by selecting a sufficiently dense grid |[T4l . |[T5l . Under 

(i) fm''') 

asl), each grid point i is associated with a unique target index at time k; that is, z G , where 

m[*^ G [1)-^] if it is occupied by one of the M targets; or, irij^ = if it is not occupied, meaning it 
is associated with a dummy target m = with strength s^^^ = 0. Apparently, {Q\^^}m=o mutually 
exclusive across m and their union spans the entire grid in the sense U^^qQ]^^ = U^^i, which leads 
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to a measurement equation [of. ([T}] 

M 

yn,k =Y. ^ (fi^^^")) X^^ + Vn,k = h^Xfc + Vn,k- (2) 

Here := [h{(f^'^),h{(f^'^), . . . ^hicf^"^)]; d*^" := ||q„ - g^Ha now denotes the known time- 
invariant distance between the nth sensor and the ith grid point; and the noise v^^k replacing i/„ ^ in 
([T]) captures the unmodeled dynamics in the aforementioned spill-over effect. Notwithstanding, thanks to 
the grid-based model, the measurements in (|2]i have become linear functions of the unknown whose 
nonzero entries reveal the grid points where target signal strengths are present at time k. 

The next step is to model the evolution of in time as the targets move across the grid. Regarding 
their movement pattern, targets obey the following assumption: 

(as2) All targets move according to identical transition probabilities where fj,''^^ := p{x^j^^ ^ 

In words, the homogeneity of targets under as2) refers to the probability that a target m moves from grid 
point i at time /c — 1 to point j at time k. 

Consider now expressing each entry of x^ as x^''^ = s^'^^ • p{x\^^ / 0), where s^^'^ = s*^™*^'') € 
|g(o) g(M)| denotes a nonnegative proportionality constant, and p{x^^^ / 0) stands for the 

probability of a target to be present on grid point j at time k. Essentially, each x): is associated with 
only one of the (M + 1) targets (including the dummy target m = 0) indexed by m\^\ and is a 
proportionality constant in the sense that it takes on (M + 1) possible values s^™) = X]jgg<'"' ^if^' 
m = 0, 1, . . . M. 

Under asl) and as2), it is shown in the Appendix that the state obeys the following recursion 

4^'^ = i:/f^-Ei, Vj€[l,G]. (3) 

Concatenating (l3]l for j = 1, . . . , G, and Q for n = 1, . . . , A^, one arrives at the grid-based model 

Xfc = FfcXfc_i + Wfc (4a) 

Yk = Hxfc + Vfc (4b) 

where the GxG state transition matrix has its (i, j)-th entry given by f^^^', the measurement matrix is 
defined as H := [hi, . . . , h„]^; likewise for the measurement noise vector := [vi^^, . . . , uat fc]; and 
is a zero-mean process noise vector with a positive-definite covariance matrix added to account for 
both asl) and the natural non-negativity constraints on x^. whose entries represent target signal strengths 
(magnitudes or power). 
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A distinct feature of model (jj]) is that the unknown is sparse \/k, since only few of its G entries 
are nonzero (in fact exactly M nonzero entries if all the M targets are located on grid points). Although 
^ describes the linear evolution of each x/^. entry under asl), using these recursions alone does not 
guarantee that the predicted or estimated x^ adheres to asl). Indeed, starting with a target at an arbitrary 
entry in xq 7^ and running (l3]l up to a large enough k, the signal strength of this target will "spill-over" 
to all entries of x^, and will possibly overlap with other targets present. Such a state transition pattern 
is expected, because uncertainty of any dynamically evolving state grows over time if no corrections 
are made based on real-time measurements. Therefore, x^ predictions based on (|4al i will be non-sparse, 
but the true state vector x/^. at any time k is sparse with only a few nonzero entries around the target 
locations. Posterior to processing the measurements, filtered and predicted renditions of x^ will remain 
sparse as well. The noise term reflects the uncertainty in the state transition model under asl). 

This sparsity attribute will prove to be instrumental for enhancing tracking performance. Also, it is 
worth noting that the state transition matrix reflects the transition behavior of target positions only, 
without revealing full information of the target movement model that may be dependent on velocity or 
other factors as well. In fact, is derived from the target movement model but does not fully reveal it, 
which differs from most existing track state models. 

Given yi-^ := {yi, . . . , y^}, the goal of this paper is to track x^ using a sparsity -aware Kalman filter 
(KF). Since x^ represents the target signal strength on the grid (TSSG), the KF-like algorithms proposed 
in Sections JII] and |IV] will be referred to as TSSG-KF trackers, while the iterated extended Kalman 
filter (lEKF) algorithms of Section |V] will be referred to as TSSG-IEKF trackers. Having available x'^p 
estimates, and recalling that x^"'^ = s*^™*^''^p(s[,^^ 7^ 0), one can estimate the constant s^"^) capturing the 
signal strength of the m-th target at time k as 

4"^^ = E,,,r) 4^'\ vfc (5) 

and the corresponding target position vector at time k as 

pf^ = (i/4'"te,.gr) - = i,...,M. (6) 

The following remark makes useful observations regarding the position estimate in 
Remark 1. A TSSG filter for tracking x^ avoids data association, because the TSSG-based state and 
measurement equations in dUl hold for any target-grid association {^^™^}m> so long as asl) and as2) 
are satisfied. On the other hand, finding the target positions via Q requires knowledge of {^^"^^}m> and 
hence calls for associating targets with TSSG entries. Solution to such an association problem will be 
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provided in Section |Vl] Nonetheless, it is worth stressing that the association errors and resultant position 
estimation errors do not affect TSSG tracking that is independent of target position estimation, similar 
to the PHD and BOF in |l22]l and IHl, respectively. 

In addition to reduced complexity, an attractive feature of the present formulation relative to e.g., f]A\ 
is that even for finite G, there is no need to assume that targets are located on grid points since Q allows 

for interpolating the target position vectors regardless, after knowing that grid point j is associated with 

(i) 

the target mj, occupying it. The next remark is useful to further appreciate this point. 
Remark 2. Given measurements yi:fc, and supposing that the number of targets M and their signal 
strengths {s^^^ , ■ ■ ■ , s^^'^^} are known, the maximum a posteriori (MAP) and minimum mean-square error 
(MMSE) optimal trackers can be derived from a hidden Markov model (HMM) filter implementing the 
following recursions derived from B ayes' rule (cf. (134) and ( [35] ) in the Appendix) 



yi:fc-l 



P ( x^^ / 



yiifc 



p{Yu\xf / 0;5(-"'))p(4-'') ^ 0|yi;fc_i) (7) 



E. („o),p(yfcl4'V0;sK'')Mxl*V0|yi;fe-i) 



where f^^^ is the transition probability as in These HMM recursions hinge on prior knowledge of the 
target-grid association {^^'"^1^=0' which need to be figured out among a total of (M + 1)^~'^G\/{G — 
M)\ possible combinations. A large G increases grid density and hence spatial resolution, at the expense 
of increasing complexity. In addition, M and {s^'")}^^^ need to be known beforehand. 

One additional remark is now in order. 
Remark 3. Although y^ in ( |4bl ) comprises scalar measurements from geographically distributed 
sensors per time k, it is possible to form y^ with samples of the continuous-time signal received at a 
single (e.g., a radar or sonar) sensor by over-sampling at a rate faster than the rate changes, so long 
as the state-space model ^ is guaranteed to be observable (and thus is ensured to be identifiable). 

m. KF FOR Tracking TSSG 

If the non-negativity constraints for x^ were absent, the optimal state estimator for dD in the MAP, 
MMSE, or least-squares (LS) error sense would be the clairvoyant linear KF. A pertinent state estimator 
is pursued here in the presence of non-negativity constraints. Suppose that the estimate ^k-i\k-i ^^d 
its error covariance matrix Pfc_i|fc_i are available from the previous time step. At time k, the KF state 
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predictor and its error covariance are obtained as 

Xfclfc-l = Fk^k-l\k-l 

(8) 

Pfc|fc-i = FfcPfc_i|fc_iFfc + Qfc. 
For the KF corrector update, consider the LS formulation of the KF; see e.g., |[27ll . The corrector 
update can be derived as a regularized LS criterion, which will also be useful to account for the sparsity 
attribute. To show this, view ^k\k-i a noisy measurement of x^. It follows readily from dD that 
Xfcifc-i = + e,t|fc-i' where ej^\j^_^i has covariance matrix P^i^.i. Stacking ^k\k-i ^^d to form an 
augmented measurement vector, yields the following linear regression model 







Ig 




^k\k~l 








Xfc + 




Yk 




H 




Vfc 



where the augmented noise vector has block diagonal covariance matrix denoted as diag(Pjt|;,_i, Rfc). 
The weighted (W)LS estimator for this linear regression problem is given by 

Xfcifc = arg min ||xfc|fc_i - Xfc||2 _i + Hy^ - Hxfc||^-i (9) 

X;^ ^0 I A: — 1 fc 

where ||x||^ := x-^Ax. In the absence of non-negativity constraints, the optimal state corrector x^i^^, can 
be found in closed form as the cost is quadratic, and likewise its error covariance can be updated as 

Pfc|fc = Pfc|fc-i - Pfc|fc-iH^(HPfc|fc_iH'^ + Rfc)~-^HPfc|fc„i. (10) 

A gradient projection algorithm will be developed in Section JV] to solve Q under non-negativity 
constraints on the state vector. However, ([TOl i will still be used bearing in mind that this update is 
approximate now. The TSSG-KF tracker implemented by (l8Tl- (fT0l) is sparsity-agnostic, as it does not 
explicitly utilize the prior knowledge that x^ is sparse. 

IV. Sparsity-aware KF Trackers 

Taking into account sparsity, this section develops sparsity-cognizant trackers. To this end, the degree 
of sparsity quantified by the number of nonzero entries of x^, namely the iQ-novm ||xfc||o, can be used to 
regularize the LS cost of the previous section. Unfortunately, similar to compressed sensing formulations 
for solving under-determined linear systems of equations [fTOl . such a regularization results in a non- 
convex optimization problem that is NP-hard to solve, and motivates relaxing the ^o-norm with its closest 
convex approximation, namely the £i-nomi. Thus, the proposed sparsity-cognizant tracker is based on 
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the State corrector minimizing the following £i -regularized WLS cost function 

Xfcifc = arg mill J(xfc) (11) 

Xfc>0 



J(xfc) := ||xfc|fc_i - Xfe||p-i + \\yk - Hxfc|L-i + 2Afc||xfc||i. 

fc I — 1 k 

The state corrector minimizing (fTTT i. together with the co variance updatq^ in dTOl ) and the prediction 
step in ([8]l, form the recursions of the sparsity-aware TSSG-KF tracker Relevant design choices and 
algorithms for minimizing ([TTI i will be elaborated in the next subsection. 

The TSSG-KF trackers in ^ and (fTTI ) involve both prediction and correction steps, which interestingly 
can be combined into a single estimation step. Considering that both Xfc_i and x^ are sparse and 
non-negative, and combining the LS terms for both the prediction and correction steps, the following 
optimization problem arises for some non-negative Xk-i and parameters: 

Xfclfc = arg min I \\^k~i\k~i - Xfc-i||p-i + ||xfc - FfcXfc„i||2 ^ 

+ llyfc — Hxfc||^-i + AA;_i||xfc_i||i + Afe||xfc||i| . 
The performance gain of this tracker was evaluated via simulations and no substantial improvement over 
the TSSG-KF tracker was observed. For this reason, focus henceforth will be placed on the TSSG-KF 
tracker in ([TTI i. 



A. Parameter selection 

The scalar parameter A^ in (ITTI ) controls the sparsity-bias tradeoff IITtI . The corrector x^i^ becomes 
increasingly sparse as A^ increases, and eventually vanishes, i.e., = 0, when A^ exceeds an upper 
bound Afc. There are two systematic means of selecting Afc. The first one popular for variable selection 
in linear regressions is cross-validation ifTTI pp. 241-249]. The second one is the so-termed absolute 
variance deviation based selection that has been advocated in the context of outlier rejection setups |fT6l . 
Both approaches require solving ([TTI i for different trial values of Afc. Even though warm starts reduce 
the computational burden considerably, this can be certainly affordable for offline solvers of a linear 
regression problem or a fixed-interval smoothing scenario, but may incur prohibitive delays for real-time 
applications. For the tracking problem at hand, the simple rule advocated is to set A^ = aXk, where 
a G (0, 1) is a fixed scaling value to avoid the trivial solution x^i^ = 0. The bound \k is derived next. 

'a more accurate covariance update will be derived in (1281 . 
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Proposition 1. The solution to (fTTI ) reduces to Xf^^k = ^ far any scalar Afc > A^, where 

Afc = \\^k[k-i^k\k~i + H'^R-fc "^yfelloo- (13) 

Proof: Since > 0, it holds that ||xfc||i = x^l, where 1 denotes the all-one vector. Therefore, J(x) in 
(fTTI) is differentiable and results in a convex problem. The necessary and sufficient optimality condition 
states that x* is an optimum point iff (y — x*)^V J(x*) > 0, Vy > 0. For x* = 0, this condition holds 
iff VJ(x*) > 0. It then follows from ^ that 

VJ(x) = 2( - - x) - H^Rfc '(yfc - Hx) + Afcl) • (14) 

Therefore, x* = is an optimal solution iff (fT3] ) holds. □ 

B. Gradient projection algorithms 

As (fTTI ) is a convex problem, convex optimization software such as SeDuMi ||26ll can be utilized to 
solve it efficiently. In addition to these solvers, low-complexity iterative methods are developed here, by 
adopting the gradient projection (GP) algorithms in Q pp. 212-217]. Note that the proposed algorithms 
can be used to obtain the sparsity-agnostic tracker from (|9} too, since the latter is obtained by minimizing 
a special case of (fTTI ) corresponding to A^ = 0. 

At each time k, the GP is initialized with x;-|;j(0) = x^i^.i at iteration / = 0. The state corrector 
iterates from / to (/ + 1) as follows 

Xfc|fe(/ + 1) = [±k\k{l) - 7V J + (15) 

where [x]~^ denotes the projection onto the non-negative orthant, 7 is the step size, and V J is as in (fT4l ). 
Here J(xfc) is differentiable because ||xfc||i = x|^l when x^ > 0. 

While (fTSl ) amounts to a Jacobi-type iteration updating all the entries at once, one can also devise 
Gauss-Seidel variants, where entries are updated one at a time fT, pp. 218-219]. This is possible because 
the non-negative orthant is a constraint set expressible as the Cartesian product of one-dimensional sets, 
allowing entry-wise updates per iteration (/ + 1) as 

x^](/ + l) = max{0, xifJ(0-7V,j(i(3](/))} (16) 

where x^'|^(/) := |x^|^-' ^\l + l),xj^^'^\l)^ has its first (j — 1) entries already updated in the (/ + l)st 
iteration. Convergence of the iterations in (fT6) to the optimum solution of (fTTI ) is guaranteed under mild 
conditions by the results in ||21 p. 219]. Specifically, J(xfc) should be non-negative and its gradient should 
be Lipschitz continuous, both of which hold for the objective in (fTTI ). 
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Proposition 2. Any limit point of the sequence generated by (116b . with arbitrary initialization x^^^, is 
an optimal solution of (|1 II) provided that the step size 7 is chosen small enough. 

In practice, only a few gradient-projection iterations are run per time step k to allow for real-time 
sparsity-aware KF tracking. 

V. Enhanced Sparsity-aware IEKF Tracking 

The proposed sparsity-aware tracker employs the KF covariance recursion in ([TOl i to update the error 
covariance of the corrector state estimate. As it does not account for the ^i-norm regularization, this 
update is approximate. In order to incorporate the prior knowledge of sparsity when updating the corrector 
covariance, this section develops an EKF-based approach, which leads to enhanced tracking performance. 

Toward this objective, the prior information on sparsity is viewed as an extra measurement /i^ = ||xfc||o, 
rather than as a regularizing term in the LS cost function. When the number of targets M is known, 
an apparent choice is to set pk = M. Accordingly, tracking will be carried out based on an augmented 
(A^ + 1) X 1 measurement vector, given by 

Yk ■■= [Yfc m ■ 

A. Viewing sparsity as an extra measurement 

The added measurement can be modeled in a general form as 

where p(xfc) is a differentiable function approximating the sparsity-inducing ^Q-norm, and denotes 
zero-mean noise with variance af,. The noise term captures both the uncertainty in approximating ||xfc||o, 
as well as the error in attaining the desired degree of sparsity. As to p(xfc), three well-known approximants 
of the ^o-norm are the ^i-norm, the logarithm, and the inverse Gaussian functions: 

(^i-norm) p(xfc) = xf'l 
(logarithm) p(xfc) = log (a:? + s) 

(inverse Gaussian) p{yik) = Y!j=i {l - exp^-^^^^^ 

where 5 and ap are tuning parameters, and only x^ > is considered. These nonlinear functions are 
plotted along with the £o-norm function for a scalar x^ in Fig. [T] It can be seen that they all have 
relatively sharp edges around the origin to approximate the iQ-novm. 
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Adding the extra measurement fik, the state space model in (UJl is augmented to 

Xfc = FfcXfc_i (17a) 

fk = h(xfc)+Vfc (17b) 

where h(xfc) := [(Hx^)^, pfc(xfc)]-^ consists of scalar measurement functions that can be nonlinear 
in general, and Vfc := [v|^, n^]^ has covariance R^. := diag(Rfc, it^). Similar to ([TTI i. the model in ([TT] ) 
leads to a nonlinear (N)LS problem 

Xfcifc = arg min Ji(xfc) (18) 
Ji(xfc) := ||xfc|fc_i - Xfc|||-i + llyfc - Hxfc||^-i + {fik - yo(xfc))^ • 

Compared with (fTTI) . ([TSl l replaces the ^i-norm of x^ with an alternative LS-enor regularization involving 
the extra measurement which accounts for the sparsity present. Because (fTSl) directly results from ([TtT i. 
the error covariance of state estimates can be updated using the KF-like recursions developed next. 

B. lEKF algorithm for nonlinear measurement models 

Since the augmented in (llVbl ) is a nonlinear function of the wanted TSSG state, the EKF approach 
is adopted here to update the error covariance along the lines of e.g., |4, Chap. 10]. Specifically, an 
iterated (I)EKF algorithm is employed, which is tantamount to applying Gauss-Newton iterations to a 
relevant NLS regression problem 10. 

The prediction step of the lEKF is similar to KF, hence ^k\k~i ^^d P^k-i follow directly from the 
state space model in (fTTI) . and coincide with dSjl. For the correction step per time k, lEKF recursions are 
initialized with Xfc|fc(0) = ^k\k-i for ^ = 0' ^^d subsequent iterations proceed as follows [28, Appendix C] 

Xfe|fc(^ + 1) = + K(/)(yfc - h(xfc|fc_i) + *(/)(xfc|fc(/) - Xfc|fc_i)) 

^ _ ^ (19) 

K(/) = Pfc|fc_i*^(/)(*(/)Pfc|fc_i*^(/) + Rfc) ' 

where := Vh(x^.|fc(/))"^ denotes the Jacobian matrix of h(-) evaluated at ^^ki^)- After the lEKF 
iterations are completed at l = L, the corrector's error covariance matrix is updated as 

Pfc|A: = Pfc|A:-l-K(L)*(L)Pfc|fc_i. (20) 

The ensuing proposition establishes the hnk between lEKF and Gauss-Newton iterations for the related 
NLS problem. 
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Proposition 3. Consider the NLS problem [cf. (fTTl ) and (fTSl )/ 



Xfcifc = arg mill \\^k\k~i - ^k\\p-i + ||yfc - h(xfc)|||-i. (21) 

Solving (1211 ) via Gauss-Newton iterations initialized with icf,^i^{0) = xi^^j^_i, amounts to the lEKF 
recursions in ( I19I ). 

Proof: The quadratic tenns in (1211 can be rewritten as 

arg min ||g(xfc)||2 



where g(xfc) 



Xfc 

-1/2 



(22) 



^k\k^ii^k\k~i-^k) 



R 



-1/2 



(yfe - h(xfc)) 



(23) 



Gauss-Newton iterations for (|22] | become 

Xfe|fc(Z + 1) = x,|fc(/) - (*(/)*'^(0)"' *(Og(xfc|fe(/)) (24) 

where ^'(Z) := Vg(xj;.|^.(/)) is the Jacobian transpose evaluated at Xjt|fc(0- Substituting g(.) from (|23] ) into 
(l24l) . and applying the matrix inversion lemma to invert the matrix in (l24l ). yields ( fT9l ) after straightforward 
algebraic manipulations. □ 
When Gauss-Newton iterations in (l24l ) are adopted in lieu of lEKF, the resulting error covariance 
matrix is a function of Vg at the last iteration L given by 

Pfc|fc = (*(L)*^(L))~\ (25) 

The sparsity-aware EKF formulation in (fTSl ) is a special case of the general NLS problem in (1211 
corresponding to h(x/fc) := [(Hxfc)^, /9fc(xfc)]^. As a result, the error covariance for the state estimate 
of (ITSl) can be derived from (|25] | as 

Pfc|fc= (^P^|i_^+H^R^iH+^Vp(ifc|fc(L))Vp(ife|,,(L))^) . (26) 

Compared with ([TOl i for the sparsity-agnostic KF, the last summand in (l26l l captures the effect of the 
sparsity-promoting penalty term on the error covariance. 

To enforce the non-negativity constraints in ([TST l. one can project each Gauss-Newton iterate in (|24] | 
onto the non-negative orthant. Unfortunately, this may not generate a convergent sequence |17] p. 215]. 
To ensure convergence, the projection should be with respect to a different distance metric than the usual 
Euclidean distance. Upon defining B(/) := (*(Z)*^(Z))~^, one implements 



ifcifc(/ + i)= [x,|fc(0 - (*(/)*^(/)) i(Og(V(0) 



B(0 



(27) 
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where denotes projection onto the non-negative orthant, which minimizes the ||.||g distance instead 
of the usual \\.\\2- If pip^k) = which is equivalent to the ^i-norm for > 0, then (fTSl ) becomes 
convex and general-purpose convex solvers such as SeDuMi can also be utilized to solve it ||26l . 

The iterative updates in (|27| | and (l26l ). along with the prediction step ([8]l, constitute the sparsity-aware 
TSSG-IEKF tracker. 

C. Enhanced sparsity-aware KF tracker 

As a final note, the sparsity-aware TSSG-KF tracker in Section |IV] can be enhanced by also casting the 
£i -regularized WLS cost in ([TTI l as an NLS cost. The ^i-norm term in (fTTT i can be equivalently expressed 
as an extra LS error term for the extra measurement = \/2A ^/x^l + Uk, where Uk is zero-mean noise 
with unit- variance. The corresponding covariance update can be derived from (l25T l as 

^k\k = (P,Y.-1 + H^I^.- + ^11^) ' • (28) 
In all, the state update in (fTTI ). together with the prediction step in ([8]) and the refined covariance update 
in (|28] ). form the recursions of the enhanced sparsity-aware TSSG-KF tracker. 

VI. Position Estimation and Track Formation 

The TSSG filters developed so far produce a dynamic TSS map of the operational environment. Such 
information is adequate to describe the targets' distribution and spatial occupancy over the sensing field 
of interest, similar in spirit to the PHD filter which portrays the targets' intensity function and the BOF 
that depicts their occupancy map. In many tracking applications however, more informative estimates 
such as target positions and trajectories are desired. This section provides TSSG-based solutions to these 
estimation tasks too. 

For the PHD approach, methods performing these extra steps have been reported using particle PHD 
filters |[T2l . II20I . II25I . or Gaussian mixture (GM)-PHD filters |24|. Target positions are typically identified 
by peak-picking the target intensity function being tracked, and the estimated target positions are treated 
as measurements for the ensuing data association and track recovery tasks. PHD filters view each particle 
or each Gaussian component involved as a target ll22l . ll30l . and employ conventional target movement 
models to describe the state transition. As a result, most of the well-known data association methods 
can be run after PHD filtering ||3], lEl Chapters 6-7]. Examples include the auction algorithm proposed 
in II20I . and the joint probabilistic data association (JPDAF) algorithm ll23l . Likewise for the BOF, the 
target movement model is employed in updating the HMM filter as well, which makes it feasible to be 
combined with a well-established data association method such as the JPDAF 1231 . 
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In contrast, the TSSG state equation only models the dynamic behavior of the TSS distribution on the 
grid, in which grid points are not treated as targets, and hence do not directly obey the conventional target 
movement model. As remarked in Section JIl only partial information about position changes is explicitly 
captured by the state transition matrix F^, while other factors such as velocity are implicit. Due to this 
major difference, conventional data association methods cannot be directly adopted as a follow-up to 
TSSG filtering. This section develops estimators of target positions and tracks for multi-target scenarios, 
based solely on the limited information regarding target transition probabilities on the grid. 

A. Target position estimation 

Given the output x^i^ of the TSSG filter, target positions can be obtained from 1^ provided that the 
subset of grid points associated with each target is known in the form of Q^i^\ 'im. 

Starting from x^i;-, one can apply appropriate clustering techniques to identify Q^jl^\ When the number 
of targets M is known, simple parametric clustering methods such as the /c-means can be used 18] pp. 424- 
429]. When Af is unknown, one can perform joint clustering and model order selection. Such algorithms 
utilize some global model order selection criteria such as Akaike's information criterion to determine 
the best number of clusters M, as well as the clusters {^|."^^}m=i themselves OTI . Other nonparametric 
clustering methods can be employed as well, without assuming or estimating the number of clusters. For 
example, hierarchical clustering techniques either aggregate or divide the data based on some proximity 
measure, while density estimation-based nonparametric approaches identify clusters and their number 
from the modes of the empirical density function of the unknowns; see e.g., |[T8l for a survey. 

Having acquired M and {Q^}^^}m=i, and based on Q, the target positions can be obtained individually 
from the TSSG estimates on the associated clusters of grid points Vi G Q\^\ as follows: 

.M^E^e^r^^ m = l,2,...,M. (29) 

B. Position-to-track association 

Suppose that there are Mj tracks from time slot 1 up to /s — 1, and p^^i has been associated with track 
t and hence alternatively expressed as p^*!^, t = 1, . . . ,Mt. The goal of track association is to assign 
the position estimates |p[.'"^| of the M targets at time k to one of the established Mj tracks. For 

I J m=l 

clarity in exposition, suppose first that M = Mt and there is no target birth or death. This assumption 
will be removed later on. Evidently, there are Ml different assignments, which must be examined to find 
the best possible association. 
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Given yi:fc-i, the first step is to establish! a track prediction model to be used for computing the 
predicted track positions {p[,*j,,_^}f=i and their error covariances. Note from (|29^ that the target position 
estimates conditioned on the TSSG are independent of the per-sensor measurements. Hence, it suffices 
to predict {p^i^_^}t solely from the TSSG vector x^,_i|;j_i. To do so, focus on track t and form a G x 1 

vector ^k-i,t that only retains the entries of ^k~i\k~i belonging to the t-th cluster of grid points in G^^l^i', 

(i) (i) (t) (i) 

that is, x)^_i I = for j G Q^Li and x^_^ ^ = otherwise, Vj. 

Given Xfc_i t at time k — 1, the predicted TSSG belonging to track t at time k becomes 

Xfc|fc-l,i = FfcXfc_i^t 

and correspondingly, the predicted track position is 

(t) _ 1^0=1 gj^fc|fc-i.t 
Pfe|fc-1- G Aj) ■ ^^^> 

^j=l •^k\k-l,t 

The normalized quantities (S^=i t) <l30l) play the role of fractional weights when the 

corresponding grid positions gj are used to estimate the track position. Viewing p|.*|j,_^ as the weighted 
average of G position-samples {gj}^=i, it is straightforward to estimate the covariance of p^*^^,_^ using 
the sample covariance, as 

- (t) _ Sj=i ^fc|l_i,t(gj - Pfc|i_i)(gi - pi|i_i)'^ 

Z^i=l •^k\k-l,t 

The process in (l30l)-(l3Tll is repeated for all target tracks t = 1,. . . ,Mt, so that the prediction estimates 
and covariances become available for all tracks. 

Now, the aim is to associate the predicted track positions {p^f,_-^}t in ( l30l ) with the target position 
estimates {p[,*^}m in (l29l ). To this end, define the decision variables a{t, m) G {0, 1} for t = 1, . . . , Mt 
and m = 1, . . . , M, where a(t, m) = 1 amounts to deciding that target m measured at p^™^ is assigned 
to track t. The pairwise-association cost can be quantified using the Mahalanobis distance between track 
t's prediction and p^'"^ as a measurement, that is 

MD(t, ,n) := (p^ - pf^)^ (Piti)"^ (4-1 " ' ^^^^ 

The following optimization problem is formulated to minimize the total association cost subject to 
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linear constraints that ensure one-to-one track-to-measurement mapping: 

M M 

mm 




(33) 



M 



such that 



^ a{t,m) 



1, Vt = l,...,Mt 



m=l 



M 



1, Vm= 1,...,M. 



t=i 



It is worth mentioning that (|33] ) is a special case of the so called assignment problem, which is a well- 
known data association algorithm ||9] pp. 342-349]. Its solution can be efficiently computed in polynomial 
time using integer programming solvers such as the Hungarian algorithm |[T9l . 

The track association problem in (|33T l can be modified to handle track birth and death scenaria ||9]. 
Toward this objective, introduce a dummy target m = and a dummy track t = 0. The one-to- 
one constraints in (1331 ) are modified as follows: each track is assigned to at most one target position 
measurement, but the dummy track can be associated with any number of targets; meanwhile, each 
position measurement is assigned to at most one track, but the dummy measurement can be assigned to 
multiple tracks; further, the dummy target cannot be associated with the dummy track. Such a modified 
association problem resembles the auction algorithm ||9], ||20l . along with the corresponding association 
costs defined in (l32l ). The computational burden of this combinatorial problem can be reduced by 
removing some unlikely association pairs in advance. Essentially, if for a track t all the association 
costs {MD(t, m)}m exceed a large threshold, then this track is considered "dead", and is associated with 
the dummy target. Similarly for a target m, if all the association costs {MD(i,m)}t are too large, then 
this target is considered "born", and is associated with the dummy track. 

Once the position-to-track association is completed, velocity estimates can be obtained too. This is 
possible by subtracting target position at time k — 1 from its position at time k and dividing by the 
sampling period. 

Finally, it is worth noting that in formulating (|33T l. only the state transition probability matrix 
is needed, regardless of the underlying target movement model. It is possible however to utilize each 
target's movement model to develop other (more effective) data association schemes, and refine the track 
estimates as well. Such association and track refinement steps will take place after every TSSG update, 
using the output of the TSSG tracker to form the position-measurements (|29] l for the ensuing parallel 
target trackers, one for each target. The results will not be fed back to the TSSG trackers, thus ensuring 
resilience of TSSG estimates to data mis-association and track estimation errors. 
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VII. Numerical Tests 

Consider a 300 x 300 square-meter surveillance region along with a 10 x 10 rectangular grid with 
equally-spaced grid points. Therefore, each grid cell is of size 30 x 30. Simulations are performed for 
both single- and multi-target scenarios. 

A. Single-target case 

A single target starts at the south-west corner of the grid at time k = 1, and moves northeast according 
to a constant velocity model 



where v denotes the target's constant velocity assumed known and given by v = (15, 15) meters per 
second; Pk-i is the previous target position; = 1 is the samphng time in seconds; and represents 
modeling noise of zero-mean and variance 0"^l2- Given this model and ignoring n^, if the target starts at 
the center of the grid cell it is currently in, then at the next time instant, it will arrive at the northeast comer 
of this grid cell conjoining the north, east, and northeast grid cells. Due to the symmetrically distributed 
noise, the target will have equal probability of falling inside each of the 4 grid cells. It is assumed that an 
is small enough so that the probability of a target moving into grid cells other than its four adjacent ones 
is negligible. The resultant movement model is as follows: a target stays on the current grid point with 
probability 1/4, and moves north, east, or northeast with probability 1/4. Whenever the target moves 
outside the boundaries of the surveillance region, tracking stops. One random realization of this movement 
model is plotted in Fig. |2] and is considered for the ensuing simulations starting with the single-target 
case. The target's signal strength is s = 10, and there are = 20 sensors distributed randomly over 
the surveillance region measuring the received TSS. The measurement noise is zero-mean Gaussian 
white with unit variance. The propagation function h{x) in ([T|l is given by h{x) = c/(c + x^) for x > 0, 
where c is chosen so that h{60) = 0.5. Apparently, h{0) = 1 and h{x) is monotonically decreasing as x 
increases. 

The proposed sparsity-agnostic and sparsity-aware TSSG-KF trackers in Sections UmiTVl are em- 
ployed to estimate the target signal strengths and position vectors over time. The position estima- 
tion accuracy is measured by the average root mean-square error (RMSE) in the form of RMSE = 



y — '^k=r \\Pk — Pfclli' where Kmax is the tracking duration and Pfc is obtained as in The 
covariance matrix of the process noise is set to = Iq in ([8]l, and 1,000 Monte Carlo runs 
over the random measurement noise are performed to compute the RMSE. To shed light on the role 



Pfc = Pfc-i + vTs + rifc 
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of the £i-norm sparsity penalty terni in ([TTI l. Fig. |3] depicts the RMSE perfomiance with respect to the 
sparsity-controUing coefficient Xk as a fraction of in ([T3] ). The sparsity-agnostic tracker corresponds 
to setting A = in ([TTI ). and is also plotted for comparison. It is seen that the sparsity-aware KF tracker 
outperforms the sparsity-agnostic one for a large range of A^ 7^ values, and A^ = 0.1 Afe appears to 
yield the lowest RMSE for this test. The optimal HMM filter exhibits the best performance, but requires 
accurate knowledge of the target signal strength. 

Fig. |4] depicts the RMSE of the sparsity-aware TSSG-IEKF tracker of ([TSl l. with = 1 and for 
different values of a^- This tracker incorporates sparsity as an extra measurement, and selects the sparsity 
model p(xfc) as the ^i-norm function. Evidently, this extra measurement is effective in promoting sparsity, 
which leads to improved performance relative to the sparsity-agnostic tracker. The noise variance of the 
sparsity measurement in (|17bl i is a design parameter chosen in accordance with the sensor measurements 
(here having unit variance). As Fig. |4] indicates, there is an optimal value of ak that attains the most 
effective tradeoff between the sensor measurements and the sparsity-induced measurement. As ak becomes 
larger, the tracker collects less information from the extra measurement, and eventually becomes sparsity- 
agnostic when (Tfc is too large. On the other hand, when Uk is too small, the tracker is predominantly 
enforcing a sparse solution without considering much the sensor measurements, which also degrades 
tracking performance. 

Both sparsity-aware TSSG trackers, the TSSG-KF tracker with Afc = 0.1 Afc and the TSSG-IEKF tracker 
with at = 2, are compared in Fig. |5] in terms of their RMSE performance versus time. The curves are 
generated using 1, 000 Monte Carlo runs. These two sparsity-aware trackers exhibit similar performance, 
both outperforming the sparsity-agnostic tracker. The clairvoyant optimal HMM filter is also tested as 
the benchmark. 

Finally, Fig. |6] demonstrates the dynamic behavior of the sparsity-aware estimator in (fTTt with A^ = 
0.9Afc. Even though the sparsity-aware TSSG-KF performs worse than sparsity-agnostic TSSG-KF for 
this value of Xk, it is chosen to demonstrate how sparsity affects the tracking process. The estimated 
TSSG state vectors are depicted over time, with a circle representing a nonzero TSS at the corresponding 
grid point. The true and estimated tracks are plotted as well. For clarity, only the projection of the target 
track on the y-direction is depicted. It is seen that the "cloud" of nonzero target signal strengths follows 
the true track. The estimated target profile is seen to be indeed spatially sparse. The size of the nonzero 
support indicates the uncertainty in target position estimates, which apparently does not grow over time, 
even when using a simple grid-induced linear KF tracker to follow the state transition pattern. 
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B. Multi-target case 

Two targets are respectively located at the south-center and west-center of the grid at time k = 1. 
They start moving according to the same movement model used for the single-target case. Fig. 17] plots 
one random realization of these target trajectories used for the ensuing multi-target test cases. Adhering 
to asl), these two trajectories do not overlap on the same grid point at the same time. The target signal 
strengths are set to be s^^^ = s^'^^ = 10. It is assumed that the trackers know the number of targets unless 
otherwise stated. There are 100 sensors deployed randomly over the surveillance region to measure the 
total received signal strengths. 

First, the position estimation method presented in Subsection IVI-AI is tested. Fig. |7] depicts the position 
estimates as circles along with the true target trajectories, for both the sparsity-agnostic TSSG-KF and the 
sparsity-aware TSSG-KF trackers with = O.lAfc. When the ^i-norm sparsity-promoting regularization 
term is not present (cf. Fig. |7]i, position estimates are rather inaccurate and some of them fall far from 
either of the two targets. In contrast, the sparsity-aware TSSG-KF in Fig. [7] results in quite accurate 
position estimates. One can clearly associate each position estimate with one of the two targets, and 
readily visualize target tracks from the position estimates. Before the position estimates are associated 
with individual targets, a pertinent performance metric quantifying estimation accuracy is the so-called 
Wasserstein distance (WD) that measures the distance between two finite sets |[T3l . Let = {p^^^jm 
denote the finite set of the true target positions at time k and Pk = {p^"^}n the set of position estimates, 
respectively. Let d{., .) stand for the Euclidean ^2-iiorm, and | • | for set cardinality. The WD between 
these two sets is defined as 

d^{Pk,Pk) = min (Ep(™)eP,Ep(.,eP.^™-^(P^'"^P^"T)'^' 
subject to Sj^^i = O^, Vn = 1, . . . , | A| 

Sn=l ^mn = Jp^ , Vm = 1 , . . . , | Pfc | . 

Fig. [8] depicts the WD for both sparsity-aware TSSG-KF and TSSG-IEKF trackers, in comparison 
with the sparsity-agnostic TSSG-KF tracker. The TSSG-IEKF tracker is implemented with fi^ = 2 and 
(Tfe = 2. The WD is evaluated by averaging over 1,000 Monte Carlo runs for each tracker. Evidently, 
both sparsity-aware designs are effective and improve the WD performance. 

The track formation algorithm of Subsection IVI-BI is investigated next for the same target realization. 
The target tracks formed using the position estimates of a single Monte Carlo run are plotted in Fig. |9] 
for the sparsity-agnostic TSSG-KF tracker. The estimated track for target 1 is not even plotted because it 
deviates too much from the true trajectory. The estimated track for target 2 shows some erratic behavior. 
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As will be discussed shortly, the unsatisfactory performance is not due to the proposed track formation 
algorithm itself; rather, it is a manifestation of inaccurate clustering that results from badly shaped TSSG 
estimates to begin with. The accuracy of the TSS map provided by the TSSG filters is essential in 
ensuring good performance of position estimates and track formation algorithms. Fig. |9] illustrates the 
track estimates obtained after processing the sparsity-aware TSSG-KF output. It can be seen that both 
targets are closely tracked. To compare these methods quantitatively, the RMSE curves for the two targets 
are plotted versus time in Fig. [TOl for 1, 000 Monte Carlo runs. It is evident that exploitation of sparsity 
markedly improves performance of the TSSG filters. In addition, sparsity-aware TSSG-KF seems to 
outperform the TSSG-IEKF for this specific setting and choice of parameters. 

To further illustrate the importance of TSSG estimation for subsequently forming position and track 
estimates. Fig. 11 la| depicts two snapshots of the TSSG heat maps after the KF prediction and correction 
steps at times k = 2 and 3. For the sparsity-agnostic TSSG-KF tracker, the correction heat map at k = 2 
seems to contain three clusters while there are only two targets. In the correction heat map at A; = 3, 
there is a single point in the lower right which is nonzero and far from both targets. This spurious point 
can have a detrimental effect during the clustering phase as it can greatly shift mean positions of the two 
clusters. These malign effects do not show up in the TSSG heat maps for the sparsity-aware TSSG-KF 
in Fig. Illbl where heat maps exhibit two compact clusters in both KF correction steps. 

Lastly, simulations for an unknown number of targets are performed on a 15 x 15 grid with the true 
and estimated target tracks plotted in Fig. [12] In this setup, targets 1 and 2 begin their movement at 
time A; = 1; at A; = 5 target 3 is bom, and at A; = 10 target 1 disappears. The sparsity-aware TSSG-KF 
is utilized in both simulations. Various clustering options are available when the number of clusters is 
unknown |[3T1 . Here a simple MATLAB routine called "silhouette" is used to determine the best number 
of natural clusters in the TSS maps. After A;-means clustering is performed, silhouette returns a value 
between —1 and 1 for every point that has participated in the clustering phase. The value that silhouette 
returns measures how well every point is explained by the cluster it belongs to, compared to other clusters. 
A value close to 1 is desirable. Therefore, silhouette values averaged over the clustered points offer a 
good measure of how well clusters explain the points which belong to them. The number of clusters 
with the largest average silhouette value is selected as the most appropriate number of clusters. It can 
be seen that the three targets are accurately tracked. However, a small erroneous track emerges close to 
target 1 for two time periods. Unfortunately, performance of the case with unknown number of targets is 
not always as accurate as shown here and more than one inaccurate track may arise. On the other hand, 
when applied to the two-target example previously considered in the absence of target births or deaths. 
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the algorithm with unknown number of targets is always successful in recovering accurate target tracks. 

VIII. Conclusions 

The problem of tracking multiple targets on a plane using the superposition of their received signal 
strengths as measurements has been investigated. A grid-based state space model was introduced to 
describe the dynamic behavior of target signal strengths. This model not only renders the nonlinear 
estimation problem linear, but also facilitates incorporation and exploitation of the grid-induced sparsity 
present. Two sparsity-aware Kalman trackers were developed to exploit this sparsity attribute: TSSG- 
KF promoting sparsity of the state estimates through £i-norm minimization, and TSSG-IEKF effecting 
sparsity by viewing it as an extra measurement. To address the challenge of updating the state estimation 
error covariances under sparsity constraints, a novel approach based on iterative extended KF and 
measurement augmentation was also developed to provide tractable and accurate covariance updates. 
Position estimation and position-to-track association issues were considered as well. The proposed trackers 
do not require knowing the number of targets or their signal strengths, and considerably reduce complexity 
when compared to the optimal hidden Markov model filter. They offer improved tracking performance 
at reduced sensing and computational cost, especially when compared to sparsity-agnostic trackers. 
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Appendix 

From the total probability argument, it holds that 

V {^^ / |, G )=Y^p {xf + 0, 41, / 0, ^ G Q^X \2 ^ 
which leads to the following equality after invoking as2) in Bayes' mien: 



P ( xf / 



,(m) 



i=l 



(34) 

^It holds trivially for the dummy target m = as well, because p{x^^^ 7^ 0|j € Q'^^) = and p{x"^^ 7^ 0, j G S'k ^) = 0. 
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Any grid point j = 1, . . . , G with a nonzero 7^ is associated with a single target index S 

[1,M] at time k, which means p(x]^' / 0,j G ' ) / for m^-'^ G [1,M]; and according to asl), 

.(i) 



p(x)^> / 0, J G ) = 0, Vm / m^-'^ or m = m)!' = 0. Invoking p{x)!' / 0) = Em=o^'(^fc 7^ 0, j G 
Q\^^), and noting that p(j G 5^^*° ^) = 1, yields 



p{x^P / 0) = p{x^p / 0, J G g^'"^"^) = p{xf ^ 0|j G g^'"^"^), Vj. (35) 

Similarly for a grid point i at time (A; — 1), there exists a target index mf^^ G [0, M] such that p{xf_^^ / 
0) = p{xf_-^ / 0,i G and ^(xi;^^ / 0, i ^ ^£1''^) = 0, Vi G [1,G]. Under asl) and as2), it 

follows from dM) and (l35]l that 

G 



0-) = V (4^'^ / 0) = (4-) ^ |j G af") ) = 5: /f V (41, / 0, z G 

i=l 



=0, Vi 



=0, Vi: m<^'li=0 

= E/f^41i> V,G[1,G]. (36) 

i=l 
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x-coordinate grid 



Fig. 1: The ^o-norm and its three approximations. Fig. 2: True target track on the grid. 
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Fig. 3: Sparsity-agnostic and sparsity-aware TSSG- 
KF trackers. 
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Fig. 4: TSSG-IEKF tracker with an extra sparsity 
measurement. 
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Time Time 

Fig. 5: Comparison of TSSG-KF and TSSG-IEKF Fig. 6: Nonzero support of estimated TSSG, true, 
trackers. and estimated tracks (y-direction only). 
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Fig. 7: True tracks and position estimates for two targets: (left) sparsity-agnostic TSSG-KF tracker, (right) 
sparsity-aware TSSG-KF tracker. Circles indicate the estimated target positions. 
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Fig. 8: WD versus time for sparsity-agnostic and sparsity-aware TSSG-KF and TSSG-IEKF trackers. 
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Fig. 9: True and estimated tracks: (left) sparsity-agnostic TSSG-KF; (right) sparsity-aware TSSG-KF; 
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Fig. 10: Tracking performance for multi-target case: (a) RMSE for target 1, (b) RMSE for target 2. 
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Fig. 11: Heat map: (a) sparsity-agnostic TSSG-KF tracker, (b) sparsity-aware TSSG-KF tracker. 
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Fig. 12: True and estimated tracks with unknown number of clusters. 
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